module Bootstrap
   use HTypes
   use sysutils
   use OptimisationTypes
   use equations
   use objfunction
   use strutils
   implicit none
   
contains


subroutine OnInitiateBootstrap(Opt)
   implicit none
   type(TOptimisation) :: Opt
   character(len = 256) :: Comment
   integer :: DataSourceFileHandle ! File handle open when reading bootstrap file
	type(TObservedData), pointer :: pObsData
	type(TDataSource), pointer :: pDataSource
   
   integer :: ioerr, iRun, jRun, NumDataSources, iDataSourceType, DataSourceType, iDataSource, iYear, iAge, &
      iObsYear, iObsAge, iRnd, handle_offset, ExtFirstYear, ExtLastYear, ExtFirstAge, ExtLastAge, &
      NumExtDataRuns, NumRuns, iExtYear, iExtAge
   double precision :: obj, SourceVal, Noise, res, v, X 
   double precision, pointer :: TempYear(:) => null()

   open(cBootstrapResultFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootResult)), &
      action = 'write', status =  'replace')
   write(cBootstrapResultFile, "(A)") '#Bootstrap result file   '   
   write(cBootstrapResultFile, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cBootstrapResultFile, "(A)") 'One iteration contains: ObjFunc, N matrix, F matrix '   
   close(cBootstrapResultFile);
   open(cBootstrapFMeanFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootFMean)), &
      action = 'write', status =  'replace')
   write(cBootstrapFMeanFile, "(A)") '#Bootstrap F mean file   '   
   write(cBootstrapFMeanFile, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cBootstrapFMeanFile, "(A)") 'One line contains: Iteration number and Mean F over years'   
   close(cBootstrapFMeanFile);
   open(cBootstrapSSBFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootSSB)), &
      action = 'write', status =  'replace')
   write(cBootstrapSSBFile, "(A)") '#Bootstrap SSB file   '   
   write(cBootstrapSSBFile, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cBootstrapSSBFile, "(A)") 'One line contains: Iteration number and SSB over years'   
   close(cBootstrapSSBFile);
   open(cBootstrapObjFuncFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootObj)), &
      action = 'write', status =  'replace')
   write(cBootstrapObjFuncFile, "(A)") '#Bootstrap Object function file   '   
   write(cBootstrapObjFuncFile, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cBootstrapObjFuncFile, "(A)") 'One line contains: Iteration number and objective function value'   
   close(cBootstrapObjFuncFile);
   
   ! Calculate Bootstrap Data Source
   ! Calculate residuals
   NumRuns = Opt%IndexInput%NumStages
   obj = objfunc(Opt, .true.)

   ! Read file input
   handle_offset = 0
   do iDataSourceType = 0, size(c_DataSource) - 1
      DataSourceType = c_DataSource(iDataSourceType)
      NumDataSources = GetNumObservedDataElms(Opt, DataSourceType)
      if (NumDataSources == 0) cycle
      do iDataSource = 0, NumDataSources - 1
         pDataSource => GetDataSource(Opt, DataSourceType, iDataSource)
         pObsData => Opt%ObservedData(DataSourceType)%Elm(iDataSource)
         if (associated(pDataSource%BootObs)) then
            deallocate(pDataSource%BootObs) 
         endif
         allocate(pDataSource%BootObs(0 : NumRuns - 1, 0 : pObsData%NumYears - 1, &
            0 : pObsData%NumAges - 1))
         if (pDataSource%BootType == btExtData) then
            open(cBootstrapFile, file = GetWorkFileName(Opt, &
               lrtrim(pDataSource%BootFile)), action = 'read', iostat = ioerr, status = 'old')
            read(cBootstrapFile, '(A)') Comment
            read(cBootstrapFile, *) ExtFirstYear, ExtLastYear
            read(cBootstrapFile, *) ExtFirstAge, ExtLastAge
            read(cBootstrapFile, *) NumExtDataRuns
            if(NumExtDataRuns > NumRuns) then
               NumExtDataRuns = NumRuns
            endif
            do iRun = 0, NumExtDataRuns - 1
               do iObsYear = 0, pObsData%NumYears - 1
                  do iObsAge = 0, pObsData%NumAges - 1
                     pDataSource%BootObs(iRun, iObsYear, iObsAge) = 0.0d0
                  enddo
               enddo
            enddo
            allocate(TempYear(ExtFirstAge : ExtLastAge))
            do iRun = 0, NumExtDataRuns - 1
               read(cBootstrapFile, '(A)') Comment
               do iExtYear = ExtFirstYear, ExtLastYear
                  read(cBootstrapFile, *) TempYear(ExtFirstAge : ExtLastAge)
                  do iExtAge = ExtFirstAge, ExtLastAge
                     iObsYear = iExtYear - pObsData%FirstYear
	                  iObsAge = iExtAge - pObsData%FirstAge
                     V =  TempYear(iExtAge)
	                  if ((0 <= iObsYear) .and. (iObsYear <= pObsData%NumYears - 1) .and. &
		                  (0 <= iObsAge) .and. (iObsAge <= pObsData%NumAges - 1)) then
                        pDataSource%BootObs(iRun, iObsYear, iObsAge) = V                           
                     endif
                  enddo
               enddo
            end do
            close(cBootstrapFile)
            deallocate(TempYear)
            if (NumExtDataRuns < NumRuns) then
               do iRun = NumExtDataRuns, NumRuns - 1
                  jRun = modulo(iRun, NumExtDataRuns)
                  do iObsYear = 0, pObsData%NumYears - 1
                     do iObsAge = 0, pObsData%NumAges - 1
                        pDataSource%BootObs(iRun, iObsYear, iObsAge) = pDataSource%BootObs(jRun, iObsYear, iObsAge)
                     enddo
                  enddo
               enddo
            endif
         else
            do iRun = 0, NumRuns - 1 
               do iObsYear = 0, pObsData%NumYears - 1
                  do iObsAge = 0, pObsData%NumAges - 1
                     iYear = iObsYear + pObsData%FirstYear - Opt%IndexInput%FirstYear 
	                  iAge = iObsAge + pObsData%FirstAge - Opt%IndexInput%FirstAge 
	                  if ((0 <= iYear) .and. (iYear <= NumYears(Opt%IndexInput) - 1) .and. &
		                  (0 <= iAge) .and. (iAge <= NumAges(Opt%IndexInput) - 1)) then
                        select case(pDataSource%BootType)
                           case (btOrgRes, btOrgPar)
                              SourceVal = pObsData%Obs(iObsYear, iObsAge)
                           case (btModRes, btModPar)
                              SourceVal = GetObsModDataValue(Opt, DataSourceType, iDataSource, iYear, iAge, .true.)
                        end select
                        select case(pDataSource%BootType)
                           case (btOrgRes, btModRes)
                              iRnd = floor(rndn(Opt%seed) * pObsData%NumPositiveUnweightedResiduals)
                              Noise =  exp(pObsData%ResArray(iRnd))
                           case (btOrgPar, btModPar)
                              X = snrn(rndn(Opt%seed))
                              Noise = exp(X * pDataSource%Sigma)
                              !write(*,*) "X, Sigma, Noise", X, pDataSource%Sigma, Noise 
                        end select
                        SourceVal = max(0.0d0, SourceVal * Noise)
                        pDataSource%BootObs(iRun, iObsYear, iObsAge)  = SourceVal
	                  end if
                  enddo
               enddo
            enddo
         endif
      enddo
   enddo
   close(cFileUnit)
end subroutine OnInitiateBootstrap


subroutine OnBootstrapStart(Opt)
   ! Copy from boostrap data to observed data
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   
   integer :: NumDataSources, iDataSourceType, DataSourceType, iDataSource
   integer :: iObsYear, iObsAge
	type(TObservedData), pointer :: pObsData
	type(TDataSource), pointer :: pDataSource
   ! Apply changes to the observed data
   do iDataSourceType = 0, size(c_DataSource) - 1
      DataSourceType = c_DataSource(iDataSourceType)
      NumDataSources = GetNumObservedDataElms(Opt, DataSourceType)
      if (NumDataSources == 0) cycle
      do iDataSource = 0, NumDataSources - 1
         pObsData => Opt%ObservedData(DataSourceType)%Elm(iDataSource)
         pDataSource => GetDataSource(Opt, DataSourceType, iDataSource)
         do iObsYear = 0, pObsData%NumYears - 1
            do iObsAge = 0, pObsData%NumAges - 1
               pObsData%Obs(iObsYear, iObsAge) = pDataSource%BootObs(Opt%Loc%Stage - 1, iObsYear, iObsAge)
            enddo
         enddo
      enddo
   enddo
   ! Reload the parameters
   if(.not.LoadParameters(Opt)) return
end subroutine OnBootstrapStart


subroutine OnBootstrapStop(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iAge
   double precision :: obj
   character(len=20) :: s
   ! Stream out bootstrap result
   open(cBootstrapResultFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootResult)), &
      action = 'write', position =  'append')
   write(cBootstrapResultFile, "(A)") 'Iteration ' // lrtrim(inttostr(Opt%Loc%Stage))
   write(s, "(F15.5)") Opt%Loc%ObjFunc
   write(cBootstrapResultFile, "(A)") lrtrim(s)
   do iAge = 0, NumAges(Opt%IndexInput) - 1
      write(cBootstrapResultFile,"(1000f10.4)") Opt%PopModel%Nay(:, iAge)
   enddo
   write(cBootstrapResultFile, "(A)") ""
   do iAge = 0, NumAges(Opt%IndexInput) - 1
      write(cBootstrapResultFile,"(1000f10.4)") Opt%PopModel%Fay(:, iAge)
   enddo
   write(cBootstrapResultFile, "(A)") ""
   close(cBootstrapResultFile);

   open(cBootstrapFMeanFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootFMean)), &
      action = 'write', position =  'append')
   write(cBootstrapFMeanFile,'(I4, 1000f7.3)') Opt%Loc%Stage, Opt%PopModel%StockSummary(:, ssFMean);
   close(cBootstrapFMeanFile);

   open(cBootstrapSSBFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootSSB)), &
      action = 'write', position =  'append')
   write(cBootstrapSSBFile,'(I4,1000f15.3)') Opt%Loc%Stage, Opt%PopModel%StockSummary(:, ssSSB)
   close(cBootstrapSSBFile)
   
   open(cBootstrapObjFuncFile, file = GetWorkFileName(Opt, Opt%IndexInput%OutFiles(roBootObj)), &
      action = 'write', position =  'append')
   write(cBootstrapObjFuncFile,'(I4, f16.6)') Opt%Loc%Stage, Opt%Loc%ObjFunc;
   close(cBootstrapObjFuncFile);
end subroutine OnBootstrapStop


end module Bootstrap